plot.ci2 <- function (x, qi = "ev", var = NULL, ..., main = NULL, sub = NULL, 
                      xlab = NULL, ylab = NULL, xlim = NULL, ylim = NULL, legcol = "gray20", 
                      col = NULL, leg = 1, legpos = NULL, ci = c(80, 95, 99.9)) 
{
  if (length(ci) < 3) {
    ci <- rep(ci, 3)
  }
  if (length(ci) > 3) {
    ci <- ci[1:3]
  }
  ci <- sort(ci)
  if (!"pooled.sim" %in% class(x)) {
    something <- list(x = x)
    class(something) <- "pooled.sim"
    attr(something, "titles") <- x$titles
    x <- something
  }
  xmatrix <- matrix(NA, nrow = length(x), ncol = length(x[[1]]$x$data))
  for (i in 1:length(x)) {
    xmatrix[i, ] <- as.matrix(x[[i]]$x$data)
  }
  if (length(x) == 1 && is.null(var)) {
    warning("Must specify the `var` parameter when plotting the confidence interval of an unvarying model. Plotting nothing.")
    return(invisible(FALSE))
  }
  if (is.character(var)) {
    if (!(var %in% names(x[[1]]$x$data))) {
      warning("Specified variable for confidence interval plot is not in estimated model.  Plotting nothing.")
      return(invisible(FALSE))
    }
  }
  if (is.null(var)) {
    each.var <- apply(xmatrix, 2, sd)
    flag <- each.var > 0
    min.var <- min(each.var[flag])
    var.seq <- 1:ncol(xmatrix)
    position <- var.seq[each.var == min.var]
  }
  else {
    if (is.numeric(var)) {
      position <- var
    }
    else if (is.character(var)) {
      position <- grep(var, names(x[[1]]$x$data))
    }
  }
  position <- min(position)
  xseq <- xmatrix[, position]
  xname <- names(x[[1]]$x$data[position])
  ev1 <- NULL
  if (qi == "pv") {
    request <- "Predicted Values: Y|X"
    if (!is.null(x[[1]]$x1)) {
      ev1 <- simulation.matrix(x, "Predicted Values: Y|X1")
    }
  }
  else if (qi == "fd") {
    request <- "First Differences: E(Y|X1) - E(Y|X)"
  }
  else {
    request <- "Expected Values: E(Y|X)"
    if (!is.null(x[[1]]$x1)) {
      ev1 <- simulation.matrix(x, "Expected Values: E(Y|X1)")
    }
  }
  ev <- simulation.matrix(x, request)
  if (is.null(ylab)) {
    ylab <- request
  }
  ci.upper <- function(x, alpha) {
    pos <- max(round((1 - (alpha/100)) * length(x)), 1)
    return(sort(x)[pos])
  }
  ci.lower <- function(x, alpha) {
    pos <- max(round((alpha/100) * length(x)), 1)
    return(sort(x)[pos])
  }
  k <- ncol(ev)
  n <- nrow(ev)
  if (is.null(col)) {
    myblue1 <- rgb(100, 149, 237, alpha = 50, maxColorValue = 255)
    myblue2 <- rgb(152, 245, 255, alpha = 50, maxColorValue = 255)
    myblue3 <- rgb(191, 239, 255, alpha = 70, maxColorValue = 255)
    myred1 <- rgb(237, 149, 100, alpha = 50, maxColorValue = 255)
    myred2 <- rgb(255, 245, 152, alpha = 50, maxColorValue = 255)
    myred3 <- rgb(255, 239, 191, alpha = 70, maxColorValue = 255)
    col <- c(myblue1, myblue2, myblue3, myred1, myred2, myred3)
  }
  else {
    if (length(col) < 6) {
      col <- rep(col, 6)[1:6]
    }
  }
  form.history <- function(k, xseq, results, ci = c(80, 95, 
                                                    99.9)) {
    history <- matrix(NA, nrow = k, ncol = 8)
    for (i in 1:k) {
      v <- c(xseq[i], median(results[, i]), ci.upper(results[, 
                                                             i], ci[1]), ci.lower(results[, i], ci[1]), ci.upper(results[, 
                                                                                                                         i], ci[2]), ci.lower(results[, i], ci[2]), ci.upper(results[, 
                                                                                                                                                                                     i], ci[3]), ci.lower(results[, i], ci[3]))
      history[i, ] <- v
    }
    if (k == 1) {
      left <- c(xseq[1] - 0.5, median(results[, 1]), ci.upper(results[, 
                                                                      1], ci[1]), ci.lower(results[, 1], ci[1]), ci.upper(results[, 
                                                                                                                                  1], ci[2]), ci.lower(results[, 1], ci[2]), ci.upper(results[, 
                                                                                                                                                                                              1], ci[3]), ci.lower(results[, 1], ci[3]))
      right <- c(xseq[1] + 0.5, median(results[, 1]), ci.upper(results[, 
                                                                       1], ci[1]), ci.lower(results[, 1], ci[1]), ci.upper(results[, 
                                                                                                                                   1], ci[2]), ci.lower(results[, 1], ci[2]), ci.upper(results[, 
                                                                                                                                                                                               1], ci[3]), ci.lower(results[, 1], ci[3]))
      v <- c(xseq[1], median(results[, 1]), ci.upper(results[, 
                                                             1], ci[1]), ci.lower(results[, 1], ci[1]), ci.upper(results[, 
                                                                                                                         1], ci[2]), ci.lower(results[, 1], ci[2]), ci.upper(results[, 
                                                                                                                                                                                     1], ci[3]), ci.lower(results[, 1], ci[3]))
      history <- rbind(left, v, right)
    }
    return(history)
  }
  history <- form.history(k, xseq, ev, ci)
  if (!is.null(ev1)) {
    history1 <- form.history(k, xseq, ev1, ci)
  }
  else {
    history1 <- NULL
  }
  if (k == 1) {
    k <- 3
  }
  all.xlim <- if (is.null(xlim)) 
    c(min(c(history[, 1], history1[, 1])), max(c(history[, 
                                                         1], history1[, 1])))
  else xlim
  all.ylim <- if (is.null(ylim)) 
    c(min(c(history[, -1], history1[, -1])), max(c(history[, 
                                                           -1], history1[, -1])))
  else ylim
  if (is.null(xlab)) 
    xlab <- paste("Range of", xname)
  if (is.null(ylab)) 
    ylab <- "Expected Values: E(Y|X)"
  par(bty = "n")
  plot(x = history[, 1], y = history[, 2], type = "l", xlim = all.xlim, 
       ylim = all.ylim, main = main, sub = sub, xlab = xlab, 
       ylab = ylab, xaxt = "n")
  polygon(c(history[, 1], history[k:1, 1]), c(history[, 7], 
                                              history[k:1, 8]), col = col[3], border = "white")
  polygon(c(history[, 1], history[k:1, 1]), c(history[, 5], 
                                              history[k:1, 6]), col = col[2], border = "gray90")
  polygon(c(history[, 1], history[k:1, 1]), c(history[, 3], 
                                              history[k:1, 4]), col = col[1], border = "gray60")
  polygon(c(history[, 1], history[k:1, 1]), c(history[, 7], 
                                              history[k:1, 8]), col = NA, border = "white")
  if (!is.null(ev1)) {
    lines(x = history1[, 1], y = history1[, 2], type = "l")
    polygon(c(history1[, 1], history1[k:1, 1]), c(history1[, 
                                                           7], history1[k:1, 8]), col = col[6], border = "white")
    polygon(c(history1[, 1], history1[k:1, 1]), c(history1[, 
                                                           5], history1[k:1, 6]), col = col[5], border = "gray90")
    polygon(c(history1[, 1], history1[k:1, 1]), c(history1[, 
                                                           3], history1[k:1, 4]), col = col[4], border = "gray60")
    polygon(c(history1[, 1], history1[k:1, 1]), c(history1[, 
                                                           7], history1[k:1, 8]), col = NA, border = "white")
  }
  if (is.null(legpos)) {
    if (leg == 1) {
      legpos <- c(0.91, 0.04, 0.2, 0.05)
    }
    else if (leg == 2) {
      legpos <- c(0.09, 0.04, 0.2, 0.05)
    }
    else if (leg == 3) {
      legpos <- c(0.09, 0.04, 0.8, 0.05)
    }
    else {
      legpos <- c(0.91, 0.04, 0.8, 0.05)
    }
  }
  lx <- min(all.xlim) + legpos[1] * (max(all.xlim) - min(all.xlim))
  hx <- min(all.xlim) + (legpos[1] + legpos[2]) * (max(all.xlim) - 
                                                     min(all.xlim))
  deltax <- (hx - lx) * 0.1
  my <- min(all.ylim) + legpos[3] * min(max(all.ylim) - min(all.ylim))
  dy <- legpos[4] * (max(all.ylim) - min(all.ylim))
  lines(c(hx + deltax, hx + 2 * deltax, hx + 2 * deltax, hx + 
            deltax), c(my + 3 * dy, my + 3 * dy, my - 3 * dy, my - 
                         3 * dy), col = legcol)
  lines(c(hx + 3 * deltax, hx + 4 * deltax, hx + 4 * deltax, 
          hx + 3 * deltax), c(my + 1 * dy, my + 1 * dy, my - 1 * 
                                dy, my - 1 * dy), col = legcol)
  lines(c(lx - deltax, lx - 2 * deltax, lx - 2 * deltax, lx - 
            deltax), c(my + 2 * dy, my + 2 * dy, my - 2 * dy, my - 
                         2 * dy), col = legcol)
  lines(c(lx - 5 * deltax, lx), c(my, my), col = "white", lwd = 3)
  lines(c(lx - 5 * deltax, lx), c(my, my), col = legcol)
  lines(c(lx, hx), c(my, my))
  polygon(c(lx, lx, hx, hx), c(my - 3 * dy, my + 3 * dy, my + 
                                 3 * dy, my - 3 * dy), col = col[3], border = "white")
  polygon(c(lx, lx, hx, hx), c(my - 2 * dy, my + 2 * dy, my + 
                                 2 * dy, my - 2 * dy), col = col[2], border = "gray90")
  polygon(c(lx, lx, hx, hx), c(my - 1 * dy, my + 1 * dy, my + 
                                 1 * dy, my - 1 * dy), col = col[1], border = "gray60")
  polygon(c(lx, lx, hx, hx), c(my - 3 * dy, my + 3 * dy, my + 
                                 3 * dy, my - 3 * dy), col = NA, border = "white")
  text(lx, my, labels = "median", pos = 2, cex = 0.5, col = legcol)
  text(lx, my + 2 * dy, labels = paste("ci", ci[2], sep = ""), 
       pos = 2, cex = 0.5, col = legcol)
  text(hx, my + 1 * dy, labels = paste("ci", ci[1], sep = ""), 
       pos = 4, cex = 0.5, col = legcol)
  text(hx, my + 3 * dy, labels = paste("ci", ci[3], sep = ""), 
       pos = 4, cex = 0.5, col = legcol)
}